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Abstract 

^  This  is  the  final  scientific  report-Jor  Grant  AFOSR- 83-0097  titled  Use  of  linear  and  nonlinear  algebraic 
methods  for  solving  large  stiff  systems  of  ODEs,  which  is  a  continuation  of  research  supported  under  Grant 
AFOSR-81-0193  titled  An  Investigation  of  the  Use  of  Iterative  Linear  Equation  Solvers  in  Codes  for  Large 
Stiff  Systems  of  ODEs.  In  this  reporVwe-sum  man  teethe  main  results  accomplished  in  this  research. 
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1.  Introduction  Matthew  j.  •  ^ 

Chief,  Technical  Information  Division 

One  of  the  major  uses  of  the  computer  in  science  and  engineering  is  the  simulation  of  physical  phenomena. 
Computer  experiments  can  be  carried  out  on  the  computer,  either  to  verify  the  theoiy  or  to  make  predictions, 
before  actual  and  costly  physical  experiments  are  attempted.  The  resulting  savings  in  cost  and  increased 
flexibility  greatly  extends  the  kinds  of  experiments  that  can  be  carried  out.  This  technique  is  heavily  used  in 
every  area  of  science  and  engineering. 

The  basis  of  such  a  simulation  is  a  mathematical  model  of  the  physical  phenomena.  Very  often,  this  takes  the 
form  of  differential  equations.  Thus  the  efficient  solution  of  these  equations  is  critical  for  the  success  of  the 
simulation  approach.  Our  research  deals  with  the  class  of  stiff  ordinary  differential  equations,  which  usually  arise 
in  modelling  physical  systems  with  many  different  time  scales.  They  are  especially  difficult  to  handle  efficiently 
on  a  computer  because  the  existence  of  wiedely  varying  time  scale  requires  high  overall  resolution  and  a 
consequent  high  cost.  Implicit  methods  are  designed  to  partially  deal  with  this  time  scale  problem  but  they  in 
turn  give  rise  to  large  systems  of  linear  equations  which  themselves  are  costly  to  solve.  The  purpose  of  this 
research  is  to  investigate  the  use  of  iterative  methods  for  solving  these  systems,  which  should  reduce  both  the 
computational  time  and  storage  requirements  over  direct  methods  of  solution.  The  end  result  is  a  computer 
program  that  can  bundle  very  large  systems  of  stiff  ODEs  which  in  *urn  allows  simulation  of  hrge  physcial 
systems. 

We  investigated  the  use  of  variants  of  the  preconditioned  conjugate  gradient  (PCG)  method  for  solving  the 
linear  equations,  and  truncated  Neuion-like  iterations  for  solving  the  nonlinear  equations.  Our  main  objective 
is  to  carry  out  a  systematic  study  of  the  methods  involved  and  to  produce  a  well-documented  and  well-tested 
computer  code  for  solving  large  systems  of  stiff  ODEs  that  incorporates  the  results  of  our  study. 

A  summary  of  the  major  results  will  be  described  in  the  next  few  sections.  A  list  of  publications  supported  by 
this  grant  is  given  in  Section  6. 

2.  Alternating  Direction  Incomplete  Factorisations 

The  use  of  a  good  preconditioning  is  often  essential  for  the  successful  application  of  the  conjugate  gradient 
method.  One  of  the  most  promising  class  of  precondition! ngs  is  the  so-called  Incomplete  Factorisations,  which 
can  be  viewed  as  factoring  the  coefficient  matrix  approximately  under  the  constraint  .that  the  factors  have  to  be 
sparse.  For  elliptic  PDEs,  quite  a  lot  of  theory  has  been  developed  for  these  preconditionings.  Most  of  these 
preconditionings  have  the  property  that,  when  applied  to  the  discretisation  of  a  self-adjoint  elliptic  opertor,  the 
number  of  iterations  required  to  reduce  some  norm  of  the  error  by  a  factor  of  c  is  C^h’^logtl/c)),  where  h  is 


the  mesh  site  used  in  the  discretization.  By  modifying  one  of  these  preconditionings,  namely,  the  Dupont- 
Kendall-Rachford  preconditioning  (DKR),  we  have  been  able  to  derive  some  error  estimates  which  suggest  that 
the  number  of  iterations  can  be  reduced  to  (^h'^ogfl/f)).  Thus  this  new  preconditioning,  which  we  call 
ADDKR,  is  asymptotically  faster  than  DKR.  Moreover,  numerical  experiments  on  some  model  elliptic  problems 
indicate  that  even  for  problems  of  moderate  size  the  new  preconditioning  is  competitive  with  DKR.  We  feel 
that  this  is  an  important  discovery,  both  theoretically  and  practically,  since  no  other  known  preconditioning o  in 
thia  close  have  been  shown  to  possess  a  better  asymptotic  rate  of  convergence  than  DKR.  This  work  has  been 
published  in  the  SIAM  Journal  of  Numerical  Analysis. 

3.  Nonlinear  Preconditioning 

One  promising  idea  that  some  people  have  noticed  in  the  application  of  Krylov  subspace  methods  for  solving 
the  linear  systems  Aw  -*  b  that  arise  in  the  inner  loop  of  a  Newton-like  iterative  method  for  solving  nonlinear 
systems  F(x)  =  0  (with  A  being  the  Jacobian  of  F  at  a  point  x)  is  the  use  of  the  directional  differencing  (  F(x  + 
tu)  •  F(x)  )  /  t  for  approximating  the  matrix-vector  product  Au  in  the  CG  algorithm.  This  avoids  explicit 
evaluations  of  the  Jacobian  and  only  requires  function  evaluations.  However,  most  of  the  preconditioning 
techniques  in  use  now  are  derived  from  the  matrix  elements  explicitly.  It  is  therefore  unclear  as  to  how  to  apply 
the  preconditionings  when  the  matrix  is  not  explicitly  available.  This  issue  dees  not  seem  to  hare  been 
addressed  in  the  literature.  We  have  obtained  some  results  in  this  direction.  We  have  derived  an  algorithm  for 
preconditioning  the  CG  iteration  with  directional  differencings  that  reduces  to  the  SSOR  preconditioning  in  the 
linear  case.  It  only  requires  evaluating  the  diagonal  elements  of  the  Jacobian  which  can  also  be  approximated 
by  function  evaluations  and  can  be  easily  incorporated  into  the  CG  iteration.  We  consider  it  an  important 
discovery  because  we  think  the  use  of  directional  differencing  will  prove  to  be  very  useful  in  the  stiff  ODE 
context.  We  have  performed  some  preliminary  numerical  experiments  with  the  algorithm  and  they  indicate  that 
the  new  algorithm  is  as  effective  as  linear  preconditioning  techniques  based  on  the  explicit  Jacobian.  This  work 
has  been  published  in  the  SIAM  Journal  of  Scientific  and  Statistical  Computing. 

4.  Basic  CG  Variant  of  LSODE 

A  program  built  around  the  package  LSODE  by  Alan  Hindmarsh  has  been  completed.  It  basically  replaces 
the  sparse  matrix  solver  in  LSODE  with  the  conjugate  gradient  solver  PCGPACK  developed  at  Yale.  The  code 
can  now  take  problems  in  the  sparse  format  (same  as  the  IA-JA  format  used  in  the  Yale  Sparse  Matrix  Package 
(YSMP)  and  in  LSODES)  and  solves  the  linear  systems  that  arise  at  each  integration  step  by  variants  of  the 
preconditioned  conjugate  gradient  method.  It  computes  the  Jacobian  exactly  by  rows  and  uses  absolute  error 


4 


\_vt  v 


control.  One  of  the  more  subtle  issues  had  been  the  stopping  criterion  for  the  conjugate  gradient  method,  which 
has  to  be  related  to  the  user  supplied  error  tolerance  for  the  ODEs.  No  major  changes  in  the  ODE  strategies 
hare  been  made  in  this  code. 

We  have  performed  a  systematic  test  of  the  performance  of  our  code  on  problems  in  the  STIFF  DETEST 
program  developed  at  University  of  Toronto  for  testing  stiff  ODE  solvers.  Our  main  objective  is  to  compare  it 
to  the  sparse  option  LSODES  in  the  LSODE  package.  The  results  confirm  the  effectiveness  of  the  use  of  CG* 
like  method  in  the  stiff  ODE  context.  For  example,  we  found  that  the  CG  code  takes  only  slightly  more  steps 
and  function  evaluations  than  LSODE  to  solve  the  problems  to  a  comparable  accuracy.  Frequently,  the  new 
code  requires  only  about  one  CG  iteration  per  Newton  iteration,  which  is  extremely  encouraging  considering  that 
no  matrix  factori rations  are  needed. 

Wc  have  performed  some  tests  with  the  SSOR-preconditioned  conjugate  gradient  method  applied  to  the  heat 
equation  in  two  dimensions.  The  numerical  results  agree  very  well  with  the  theoretical  convergence  estimates. 
When  h  is  small,  only  one  or  two  CG  iterations  are  needed.  As  h  grows,  more  iterations  are  needed.  As  the 
mesh  gets  denser,  the  number  of  iterations  increases  at  the  correct  aysmptotic  rate.  Generally,  only  very  few 
(relative  to  the  sire  of  the  linear  systems)  CG  iterations  are  needed  at  each  integration  step  and  the  use  of 
preconditioning  seems  to  improve  the  efficiency.  These  re salt*  show  that  the  use  of  the  conjugate  gradient 
method  is  very  promising,  as  compared  to  direct  sparse  solvers.  A  technical  report  has  been  completed  and 
submitted  for  publication. 

5.  Comparison  of  Krylov  sub  space  methods 

A  host  of  algorithms  belonging  to  the  same  class  of  Krylov  subspace  methods  have  been  recently  at  the  focus 

of  several  researchers  dealing  with  the  solution  of  large  sparse  no nsym metric  linear  systems.  Among  them  we 

mention  the  subclass  of  the  ORTHOMlN(kXVinsome),  GCR(k)  algorithms,  (Eisenstat,  Elman  and  Schultt), 

ORTHORES(k),  ORTHODIR(k)  developed  by  Young  and  Jea,  IOM(k)  (Saad),  Axelsson’s  method,  Lanctos’ 

method  and  others.  Although  the  literature  is  rich  in  methods,  few  comparisons  of  the  performances  have  been 

« 

proposed. 

la  a  recent  paper  by  Saad  [3]  dealing  with  the  Krylov  subspace  methods,  some  comparisons,  both  theoretical 
and  practical,  have  been  established.  In  particular  it  was  shown  that  the  ORTHO$ES(k)  method  of  Jea  and 
Young  is  in  fact  equivalent  to  the  10M(k)  method,  i.e.  that  the  iterates  produced  by  both  algorithms  are 
identical.  On  the  other  hand  both  numerical  stability  and  computational  cost  are  in  favor  of  the  IOM(k) 
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algorithms.  A  number  of  comparisons  shown  in  Howard  Elman's  recent  PhD  thesis2  favor  the  Chebyshev  based 
algorithms,  in  general,  against  the  Conjugate  gradient  type  methods.  A  common  argument  raised  against  the 
Chebyshev  algorithm  is  that  it  requires  the  knowledge  of  some  of  the  eigenvalues  of  A.  While  there  is  no  doubt 
that  the  process  of  estimating  the  eigenvalues  and  obtaining  the  best  ellipse  containing  the  spectrum  of  A  may 
seem  somewhat  cumbersome,  there  are  a  few  important  points  that  can  make  the  process  highly  competitive: 

1.  Once  the  optimal  parameters  are  known  the  method  is  extremely  fast; 

2.  There  are  some  hybrid  methods  under  investigations  which  combine  the  expensive  method  (e.g. 
GCR(k))  and  the  inexpensive  Chebyshev  algorithm  in  an  efficient  way; 

3.  In  the  context  of  time  dependent  problems,  the  problem  of  estimating  the  parameters  is  eased  by 
the  fact  they  are  likely  to  evolve  smoothly  thus  allowing  to  use  a  set  of  parameters  for  a  large 
number  of  steps. 

In  solving  stiff  ODE’s,  the  choice  of  a  method  is  rendered  difficult  because  one  has  to  select  a  method  which  is 

sufficiently  robust  to  handle  problems  of  different  types  and  nonetheless  able  to  take  advantage  of  the  particular 

context  of  ODE’s.  Many  of  the  methods  described  above  may  break  down  in  the  same  way  that  the  conjugate 

* 

gradient  method  breaks  down  for  indefinite  systems.  Some  theoretical  results  show  that  for  truely  stiff 
problems,  almost  all  the  eigenvalues  will  lie  in  the  right  half  of  the  complex  plane.  However,  it  is  still  important 
to  devise  algorithms  which  are  more  robust.  We  have  developed  a  stable  version  of  the  Krylov  subspace  method 
for  indefinite  and  nonsymmetric  problems  has  been  developed.  It  is  similar  but  more  economical  than  the 
SYMMLQ  algorithm  for  indefinite  problems  and  the  same  algorithm  can  handle  nonsymmetric  problems  as  well. 
This  work  has  been  published  in  the  Siam  Journal  of  Scientific  and  Statistical  Computing. 
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